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1.0  INTRODUCTION 

The  interaction  of  an  oblique  shock  wave  with  a  boundary  layer 
is  a  phenomenon  of  considerable  interest  and  frequent  occurrence  in 
supersonic  and  hypersonic  flows.  The  effects  of  shock-wave  impinge¬ 
ment  from  externally  carried  aircraft  stores  on  the  aircraft  and  the 
effects  of  shock  impingement  from  the  aircraft  on  external  stores 
have  long  been  of  interest  to  the  military.  Vehicles  such  as  the 
TITAN  I1IC  and  the  Space  Shuttle  have  geometries  that  make  considera¬ 
tion  of  boundary-layer/ shock-wave  interaction  a  necessity.  Numerous 
experimental  studies  have  been  undertaken  to  study,  correlate,  and 
explain  the  phenomenological  aspects  of  the  interaction  problem. 

Concurrent  with  the  interest  shown  by  experimentalists  in  the 
boundary-layer /shock-wave  interaction  problem,  much  effort  has  been 
expended  to  develop  theoretical  techniques  capable  of  accurately 
predicting  the  salient  features  of  the  problem. 

The  complexity  of  the  interaction  between  a  shock  wave  and  a 
boundary  layer  gives  rise  to  phenomena  not  characteristic  of  either 
a  shock  wave  or  a  boundary  layer.  Figure  1  schematically  illustrates 
the  physics  of  a  typical  boundary- layer /shock-wave  interaction.  The 
boundary-layer  equations  are  parabolic  and  hence  can  be  integrated 
(at  least  to  a  point  near  separation)  in  a  step-by-step  downstream 
fashion  once  an  impressed  pressure  gradient  is  specified.  The 
boundary-layer/ shock-wave  interaction  problem  is  not  parabolic  since 
the  impressed  pressure  gradient  is  determined  in  part  by  the  response 
of  the  boundary  layer  to  the  shock  wave.  Hence  no  a  priori  computation  of 
the  pressure  is  possible.  Moreover,  since  separation  and  reattachment  is 
a  possibility,  conventional  bound ary- layer  methods  cannot  be  used  because 
a  square-root-type  singularity  exists  at  separation  for  the  boundary-layer 
equations.  Thus,  any  solution  technique  must  appeal  to  a  system  of 
governing  equations  more  fundamental  than  either  the  Euler  equations  for 
inviscid  flow  or  the  boundary-layer  equations.  The  more  general  Navier- 
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Distance  from  Plate  Leading  Edge 


Figure  1.  Schematic  of  boundary  -layer/shock— wave  interaction. 
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Stokes  equations,  from  which  the  Euler  equations  and  the  boundary-layer 
equations  are  derived,  are  fundamental  to'the  subject  of  viscous  fluid 
flow  and  are  valid  throughout  the  entire  flow  field.  The  Navier-Stokes 
equations  will  yield  valid  results  at  separation  and  reattachment,  within 
the  separated  recirculation  region,  across  incident  and  reflected  shock 
waves,  throughout  expansion  fans,  and  for  any  combinative  influences  of  the 
aforementioned.  The  equations  require  numerical  solution  in  either  a 
spatially  elliptic  or  a  temporally  hyperbolic  domain. 

The  present  report  compares  numerical  results  from  an  explicit 
time-dependent  compressible  Navier-Stokes  analysis  developed  by 
MacCormack  (Ref.  1)  with  experimental  data  for  hypersonic  laminar 
boundary-layer /shock-wave  interaction  on  a  flat  plate  under  AEDC  von 
K£rm5n  Gas  Dynamics  Facility  (VKF) ,  Hypersonic  Wind  Tunnel  (B) ,  condi¬ 
tions  as  well  as  NASA  Langley  Research  Center  (LRC)  Mach  8  Variable 
Density  Tunnel  conditions.  It  is  shown  that  numerical  solutions  of 
the  time-dependent  compressible  Navier-Stokes  equations  yield  reasonable 
results  when  applied  to  hypersonic  laminar  bound ary- layer /shock-wave 
interactions. 


2.0  ANALYTICAL  ANALYSIS 

2.1  TIME-DEPENDENT  METHOD  OF  SOLUTION  OVERVIEW 

The  time-dependent  method  starts  with  a  complete  specification  of 
the  flow  field  and  then  uses  the  governing  equations  of  motion  to 
advance  the  flow  field  temporally  until  a  steady  state  is  reached.  Thus, 
the  flow  field  evolves  numerically  in  a  process  analogous  with  physical 
reality.  ‘  The  initial  specification  can  be  just  a  uniform  flow  with 
appropriate  boundary  conditions  (see  Roache  (Ref.  2)  and  Richtmyer 
and  Morton  (Ref.  3)  for  general  expositions  on  the  time-dependent  method). 
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The  physical  basis  of  the  time- depen dent  approach  as  well  as  the 
advantages  accrued  by  application  of  the  time-dependent  method  to 
solve  the  compressible  Navier-Stokes  equations  were  delineated  by 
Crocco  in  1965  {Ref.  4).  Kurzrock  and  Mates  (Ref.  5)  in  1966  used 
the  time-dependent  approach  to  study  analytically  the  flow  in  a  shock 
tube  and  hypersonic  flow  over  a  sharp  flat  plate.  Skoglund  and  Gay 
(Ref.  6)  applied  the  time-dependent  method  to  the  computation  of 
laminar  boundary-layer /shock-wave  interactions  in  1969. 

2.2  THE  EXPLICIT  TIME-DEPENDENT  METHOD  OF  MACCORMACK 

The  so-called  method  of  MacCormack,  since  its  introduction  in 
1969,  has  become  one  of  the  most  widely  used  explicit  second-order 
accurate  methods  for  numerical  solution  of  hyperbolic  partial  differential 
equations.  The  algorithm  was  first  introduced  by  MacCormack  in  Ref.  7 
and  subsequently  modified  and  extended  by  Refs.  8  through  14  as  well  as 
Ref.  1.  It  has  been  applied  to  obtain  solutions  of  the  time-dependent 
compressible  Navier-Stokes  equations  by  Baldwin  and  MacCormack  (Refs. 

9  and  10),  MacCormack  and  Baldwin  (Ref,  11),  and  Deiwert  (Refs.  15  and 
16)  among  others. 


The  two-dimensional  time-dependent  Navier-Stokes  equations, 
neglecting  body  forces  and  heat  generation,  can  be  written  in  conservation 
form  as: 


where 


du  <?F  dG 

dt  <3x  t?y 


(D 


f  P 


pv 


e 


(2) 
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with 


pu 

pu2  +  ax 
Puv  +  rxy 

(e  +  <rx)u  +  ryxv  _  kdT/dx 

pv 

PUV  +  rvx 

PV2  +  £Ty 

(e  +  o  )v  +  t  u  -  k  (9T/(?v 
/ 


(see  the  Nomenclature  for  terminology) • 


(3) 

(4) 

(5) 

(6) 
(7) 


The  procedure  used  to  advance  the  dependent  variables  (p,  pu, 
pv,  e)  from  a  time  t  to  a  time  t  +  At  at  the  interior  points  will  be 
examined  first.  The  procedures  for  the  boundary  points  are  different 
and  will  be  reviewed  subsequently.  The  MacCormack  method  is  of  the 
predictor-corrector  type  and  can  be  utilized  in  such  a  manner  that  a 
single  predictor  and  a  single  corrector  application  will  advance  the 
dependent  variables  in  time  by  an  amount  At.  However,  if  the  concept 
of  splitting  is  employed,  simplicity  and  computational  efficiency 
result  (Ref.  8).  Basically  the  concept  of  splitting  involves  a 
predictor-corrector  pass  driven  by  gradients  in  the  x-direction  and  a 
separate  predictor-corrector  pass  driven  by  gradients  in  the  y-direction. 
Thus  four  sweeps,  two  predictor  and  two  corrector,  are  required.  This 
may  be  written  as: 
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and 


y  n+  1/2 


) 


y  n+ 1 / 2 


n+  1/2 

i.j+1 


*.y,/2) 


(8) 

(9) 


U^1 


u"1/2  -  (f;)(F"T'/2 


p  n+l/2\ 

H.J  / 


(10) 


i  •  n-i- 1 
Li.J 


|(u"j1/2  +  u  i"!1 


p  n+  1 


) 


(ID 


Denoting  by  Lx  the  operation  performed  by  Eqs.  (8)  and  (9)  and 
by  Ly  the  operation  performed  by  Eqs.  (10)  and  (11),  the  sequence 
becomes 

Un+1  =  LxLyU”  (12) 


This  operation  (Eq.  (12)) 


Un+1  =  Ly(^)  Lx(|i)  Lx  ^  Ly(^)lT 


(13) 


retains  second-order  accuracy.  The  stability  criterion  (the  Courant- 
Friedrichs-Lewy  condition)  for  the  y  sweep  is: 

a. 


and  for  the  x-sweep 


At  = 


(15) 
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The  minimum  of  Eqs.  (14)  and  (15)  for  the  entire  grid  represents  the 
maximum  At  for  which  computational  stability  is  ensured.  The  sequence 
of  operations  defined  by  Eq.  (13)  is  used  to  advance  temporally  the 
interior  points  of  the  computational  grid.  The  boundary  condition 
treatment  and  the  initial  condition  specifications  are  needed  to 
complete  the  method. 


The  numerical  solution  of  flows  containing  strong  shock  waves  is 
often  hampered  by  numerical  oscillations  which  can  eventually  cause 
program  failure.  A  fourth-order  damping  term,  effective  only  in 
regions  of  large  pressure  gradients,  has  been  used  to  reduce  the 
numerical  oscillations  (Refs.  9,  10,  and  11).  Essentially,  this  technique 
adds  an  additional  "viscosity"  proportional  to  the  second  derivative  of 
pressure  to  each  of  the  steps  represented  by  Eqs.  (8)  through  (11). 

The  addition  in  regions  of  low-pressure  gradients  is  negligible  and 
is  of  importance  only  where  pressure  gradients  or  pressure  oscillations 
are  large.  By  using  the  arrow  symbol  to  denote  replacement,  the  damping 
terms  can  be  included  in  Eqs.  (8)  through  (11)  by 


G"j 


pn+1/ 2 

>.j 


G?tj  +  AGfj 


pn-t- 1  /  2  .\pn-i- 


(16) 

(17) 


where 


ac;  j  .  I  (i  v|  t  «>, ,  Mf'V.M-1  (u,  1+,  -  v, .) 

*,)  2  M  Pi>j+i  +  2  P;>j  +•  Pi._l  ’’J+1  I>] 


08) 


AF"V/2  -  i(|«|  +  c),,  ■l4l'i  '  l"1’1  " 


.  .  _  u.  )  (19) 

J  Pi+l,j  +  2Pi,j  +  Pi-M  ’’  M 


The  quantities  .  ,  and  F*?+^j  are  treated  in  the  same  manner  as  indicated 
1,3-1  i-l.j 

in  Eqs.  (16)  through  (19)  as  are  all  the  barred  quantities;  i.e.,  to  each 
F  or  G  term  in  Eqs.  (8)  through  (11)  is  added  the  term  analogus  to  Eq.  (18) 
or  (19) . 
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2.3  INITIAL  AND  BOUNDARY  CONDITIONS 

Initially,  the  flow  field  must  be  specified  completely  in  the 
region  under  consideration,  the  computational  plane.  The  computational 
plane  for  the  boundary-layer/shock-wave  interaction  is  shown  schemat¬ 
ically  in  Fig.  2.  All  of  the  flow  field  except  the  upper  boundary 
(line  AB  in  Fig.  2)  need  be  specified  only  as  a  uniform  flow.  The 
upper  boundary  is  specified  in  such  a  manner  that  the  incident  shock 
impacts  the  surface  at  the  desired  location.  By  specifying  uniform  flow 
along  line  AE  and  Rankine-Hugoniot  flow  conditions  along  line  EB,  such 
can  be  accomplished.  The  boundary  layer  will  form  normally  on  the 
surface,  the  shock  will  spread  downward,  and  the  interaction  will  evolve 
in  the  correct  manner. 


Figure  2.  Computational  mesh  system  for  boundary— layer/shock— wave  interaction. 


As  the  solution  advances,  only  the  boundary  conditions  along  AB, 
AC,  and  CD  must  be  specified.  Along  AEB,  the  same  conditions  initially 
specified  must  be  used.  Segment  BD  cannot  be  defined  as  this  would 
overspecify  the  problem.  A  gradient  condition  such  as  3/3x  =  0  along 
BD  is  used  since  a  large  portion  of  the  flow  field  is  supersonic  and 
will  thus  not  propogate  errors  upstream.  This  specification  of  region 
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BD  does  not  introduce  spurious  information  into  the  remainder  of  the 
flow  field.  Region  AC  can  be  taken  as  a  uniform  flow  field  if  the 
leading  edge  of  the  plate  is  to  be  considered  or  input  as  previously 
determined  profiles  if  the  area  of  interest  is  too  large  for  normal 
processing.  Both  input  profiles  and  input  uniform  flow  have  been  used 
with  good  results.  The  wall  boundary  conditions  are  defined  along  CD 
in  such  a  manner  that  flat  plate  results  are  obtained  midway  between 
the  first  two  grid  rows;  i.e.. 


ui,l  =  ~ui,2 

(20) 

vi,l  =  ~vi,2 

(21) 

1.1  -  CV(2TW  -  e;>2/C 

(22) 

.  /dp\ 

Pi,l  =  Pi,2  -  A%A,2 

(23) 

pi,l 

Pi’1  "  <y  -  u-u 

(24) 

2.4  COMPUTATIONAL  GRID 

Some  consideration  on  the  nature  of  the  expected  solution  is 
relevant  at  this  point.  The  flow  field  may  be  viewed  as  being  composed 
of  an  outer  region  basically  inviscid  in  nature  and  an  inner  viscous- 
dominated  region.  The  outer  region  will  change  only  gradually  as  the 
Interaction  evolves;  and,  therefore,  will  need  to  be  computed  less 
frequently.  Since  the  maximum  allowable  time  step  in  an  explicit 
time-dependent  method  is  proportional  to  the  grid  spacing  (as  indicated 
by  Eqs.  (14)  and  (15)),  the  finer  the  grid  the  smaller  the  allowable 
time  step  At.  If  a  fine  grid  is  used  in  the  viscous-dominated  region 
and  a  coarse  grid  is  used  in  the  outer  inviscid-like  region,  com¬ 
putations  will  need  to  be  performed  in  the  coarse-grid  region  only  once 
for  every  M  times  in  the  fine-grid  region.  M  is  the  smallest  integer 
plus  one  in  the  quotient  of  the  smallest  At  in  the  coarse  grid  and  the 
smallest  At  in  the  fine  grid,  i.e., 

M  =  MOD(Atc,  Atf)  +  1  /25) 
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Figure  2  illustrates  the  computational  grid  used,  showing  the  fine  and 
coarse  regions.  Thus,  most  of  the  computational  effort  will  be  expended 
in  the  viscous-dominated  region.  Caution  must  be  exercised  at  the 
boundary  between  the  fine  and  coarse  grids  as  accumulation  or  origina¬ 
tion  of  conserved  quantities  can  introduce  computational  anomalies  to 
the  solution. 

The  basic  digital  computer  code  follows  the  University  of  Tennessee 
Space  Institute  Short  Course  notes  as  given  by  MacCormack  (Ref.  17). 

All  numerical  results  reported  herein  were  generated  using  an  IBM 
370/165  digital  computer.  The  computer  program  was  executed  using 
single  precision  arithmetic  and  required  372,000  bytes  of  core  for  a 
183  by  35  grid  array  using  the  IBM  FORTRAN  IV  H-LEVEL  21.7  Compiler. 
Typically  about  100  sec  of  CPU  time  were  required  for  one  complete  pass 
through  the  program  (advancement  of  the  field  from  time  t  to  time  t  + 

At)  for  a  183  by  35  array.  The  code  was  modified  to  have  restart 
capability  and  to  permit  various  initial-condition  options  to  be  exercised. 

3.0  RESULTS 

The  von  Karman  Facility  at  AEDC  is  concerned  primarily  with  high 
supersonic  and  hypersonic  flows.  Since  only  boundary- layer /shock- wave 
interactions  that  are  completely  laminar  are  examined  in  this  report, 
relatively  low-incident  shock-wave  angles  are  required  at  the  hypersonic 
Mach  numbers  of  interest.  Large-incident  shock-wave  angles  at  hyper¬ 
sonic  Mach  numbers  and  the  corresponding  large  pressure  increases  across 
the  incident/reflected  shock-wave  systems  preclude  the  maintenance  of 
laminar  flows  throughout  the  interaction  region.  Laminar  interactions 
were  considered  since  it  was  deemed  desirable  to  validate  the  explicit 
time-dependent  method  for  numerically  solving  the  Navier-Stokes  equations 
prior  to  any  considerations  of  the  complicating  effects  of  turbulent 
flow.  Attempts  at  turbulent  boundary-layer/ shock-wave  interaction 
solutions  by  Baldwin  and  MacCormack  (Refs.  9  and  10)  and  by  Horstman, 
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et  al.  (Ref.  18),  indicate  that  much  development  is  needed  in  turbu¬ 
lence  modeling  before  transitional  and  turbulent  processes  can  be 
adequately  treated. 

3.1  BOUNDARY-LAYER/SHOCK-WAVE  INTERACTIONS 

Prior  to  discussing  numerical  solutions,  it  is  appropriate  to 
examine  typical  wind  tunnel  models  and  inviscid  solutions  for  shock 
reflection.  Figure  3  {taken  from  the  study  reported  in  Ref.  19)  shows 
typical  VKF  Tunnel  B  apparatus  used  to  investigate  boundary-layer/shock- 
wave  interactions.  The  upper  wedge  (generator)  produces  an  oblique 
shock  wave  that  impinges  on  the  lower  wedge  (receiver).  Impingement 
location  as  well  as  shock  strength  can  be  easily  controlled  using  such 
an  arrangement.  Pressure  or  heat-transfer  data  are  taken  from  the 
receiver  plate.  Because  of  viscous-induced  effects  near  the  leading 
edge  of  the  generator,  the  oblique  shock  angle  produced  is  not  the 
"wedge"  shock  angle  that  the  nominal  generator  angle  of  attack  would 
produce.  The  receiver  plate  leading  edge  also  exhibits  viscous-induced 
effects,  one  obvious  result  being  the  leading-edge  shock.  The 
aforementioned  effects  are  such  that,  for  nominal  given  angles,  the 
inviscid  pressure  ratios  across  the  incident/reflected  shock  waves 
ate  not  obtained. 

Numerical  solutions  corresponding  to  particular  experimental 

cases  have  been  generated  by  using  the  experimental  ratio  p^/p.|  (see 

Sketch  1)  and  the  free-stream  Mach  number  (M  )  to  define  the  inviscid 

00 

shock  angle  (6).  This  allowed  the  correct  pressure  ratio  to  be  used 
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Figure  3.  Schlieren  photograph  of  shock  generator  and 
receiver  plate  in  AEDC  Tunnel  B  (taken  from 
the  study  reported  in  Ref.  19). 
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without  expending  additional  computational  effort  to  define  viscous- 
induced  effects  on  the  generator.  Figure  4  was  computed  using  inviscid- 
f low  wedge  results  and  allows  the  rapid  estimate  of  the  appropriate 
incident  shock  angle  (6)  for  a  given  free-stream  Mach  number  and  pressure 
ratio.  Additionally,  the  minimum  and  maximum  incident  shock  angle  for 
a  normal  reflection  is  given.  It  is  interesting  to  note  that,  for  a 
diatomic  gas  (y  =  1.4)  and  free-stream  Mach  number  greater  than  about  3, 
a  regular  reflection  is  not  admissible  for  a  shock  angle  greater  than 
40  deg.  A  more  detailed  examination  was  given  by  Zumwalt  (Ref.  20). 

Figure  5  (taken  from  the  study  reported  in  Ref.  19)  is  a  schlieren 
photograph  of  a  typical  laminar  boundary-layer/shock-wave  interaction. 

The  incident/ref lected  shock-wave  system  can  be  seen.  The  dashed 
lines  indicate  the  nominal  region  used  in  the  computation.  However, 
for  a  region  which  does  not  contain  the  leading  edge  of  the  plate, 
suitable  profiles  in  density,  velocity,  energy,  and  pressure  must  be 
available  for  use  as  initial  and  boundary  conditions  for  the  upstream 
(left-hand  side  in  Fig,  5)  computational  boundary.  The  required 
profiles  could  have  been  generated  by  a  typical  parabolic  boundary- 
layer  program,  but  the  viscous- induced  leading  edge  effects  prominent 
in  hypersonic  flow  (see  Fig.  3)  would  not  have  been  present.  Since 
a  compressible  Navier-Stokes  solution  is  capable  of  generating  such 
effects,  the  computer  code  without  an  incident  shock  wave  possessed  the 
capability  of  calculating  laminar  flat-plate  profiles  including  viscous- 
induced  leading— edge  effects.  Thus,  whenever  it  was  desired  to  start 
an  interaction  computation  without  including  the  plate  leading  edge, 
a  suitable  set  of  initial  profiles  was  obtained  by  generating  a 
flat-plate  Navier-Stokes  solution  including  the  leading  edge  and 
retaining  the  required  numerical  information  on  a  data  storage  device. 
Ealdwin  and  MacCormack  (Ref,  9)  used  this  technique  in  their  turbulent 
boundary-layer/shock-wave  calculations  as  did  Carter  (Ref.  21)  in 
his  laminar  calculations  on  a  supersonic  ramp. 
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Figure  4.  In  viscid  pressure  ratios  for  incident-reflected 
shock  wave. 
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Figure  5.  Schlieren  photograph  of  a  boundary— layer /shock— wave 
interaction  in  AEDC  Tunnel  B  (taken  from  the  study 
reported  in  Ref.  19). 
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3.2  NUMERICAL  RESULTS 


As  a  check  on  the  capability  of  the  code,  a  low  Reynolds  number  flat- 
plate  solution  was  generated  using  the  conditions  of  Carter  (Ref.  21). 

The  resulting  pressure  distribution  is  shown  in  Fig.  6  along  with  the 
time-dependent  compressible  Navier-Stokes  numerical  solution  of  Carter. 
Both  solutions  exhibit  essentially  the  same  behavior  except  very  close 
to  the  leading  edge  where  the  explicit  MacCormack  scheme  is  smoother 
than  the  Brailovskaya  scheme  used  by  Carter. 


Figure  6.  Flat  plate  pressure  distributions  computed  by  the  MacCormack  Algorithm 
and  the  Brailovskaya  Algorithm. 

Depending  on  the  state  of  the  boundary  layer  and  the  magnitude  of 
the  pressure  jump  across  the  incident/reflected  shock-wave  system,  flow 
separation  and  reattachraent  with  a  recirculation  region  may  occur.  This 
condition  is  schematically  illustrated  in  Fig.  1.  The  presence  of  a 
separation  region  with  recirculation  and  reattachment  poses  a  more 
severe  test  of  the  code's  capability  than  an  unseparated  case.  Thus, 
to  avoid  possible  complications  resulting  from  separated  regions,  the 
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first  results  examined  will  be  for  non-separated  hypersonic  flows. 

All  the  numerical  results  reported  herein  were  obtained  using  a  constant 
wall  temperature  for  the  receiver  plate.  Data  which  would  allow  a  more 
exact  specification  of  wall  temperature  were  not  generally  available. 

All  of  the  VKF  Tunnel  B  data  examined  in  the  present  study  had 
the  nominal  shock  interaction  point  one  foot  from  the  leading  edge  of 
the  plate.  The  parameters  used  in  the  numerical  solution  are  given 
in  Table  1 ,  which  contains  tabular  information  for  all  the  cases 
examined.  The  first  laminar  boundary-layer/ shock-wave  interaction 
presented  (number  2  of  Table  1)  was  generated  using  a  set  of  profiles 
(density,  velocity,  internal  energy)  computed  via  the  time-dependent  com¬ 
pressible  Navier-Stokes  code  with  no  impinging  shock  wave.  Figure  7  is  a 
partial  reproduction  of  the  input  velocity  ratio  profiles  (u/U^)  generated 
and  is  typical  of  the  profiles.  Figure  8  shows  the  computed  streamline 
shapes  for  the  leading  edge  of  the  flat  plate,  and  Fig.  9  shows  the  shape 
of  the  leading-edge  shock.  Shock-wave  locations  in  the  "shock-capturing” 
MacCormack  code  were  estimated  by  computing  the  first  and  second  derivatives 

of  pressure  with  respect  to  longitudinal  distance  and  locating  regions  where 

2  2 

dp/dx  is  a  maximum  and  d  p/dx  is  large.  This  is  similar  to  the  philosophy 
of  Grossman  and  Moretti  (Ref.  22)  in  locating  incipient  shock  waves  in 
inviscid  time-dependent  calculations. 

By  using  the  aforementioned  input  profiles  and  a  shock-wave  angle 
of  8.6  deg,  a  boundary-layer/shock-wave  interaction  was  generated. 

This  corresponded  to  conditions  for  which  VKF  Tunnel  B  data  were 
available  (Ref.  19).  Figure  10  presents  the  results  of  the  computed 
interactions  as  well  as  the  VKF  Tunnel  B  data  for  the  corresponding 
tunnel  conditions.  Agreement  is  satisfactory  particularly  near  the 
peak  pressure  location  where  the  correct  shape  is  generated  and 
in  the  mid-range  where  the  correct  pressure  gradient  is  obtained.  No 
evidence  of  separation  was  obvious  from  the  data,  and  no  separated 
region  was  predicted  by  the  program.  Because  of  the  rather  modest 
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pressure  rise  across  the  incident/reflected  shock-wave  system  and 
because  of  the  lack  of  separation,  it  is  probable  that  this  is  a 
completely  laminar  interaction.  Overall  agreement  with  the  laminar 
time-dependent  Navier-Stokes  solution  tends  to  confirm  this.  Computer 
time  on  an  IBM  370/165  was  2  hr  for  the  flat-plate  solution  from  which 
the  initial  profiles  were  secured  and  3  hr  for  the  interaction. 


Table  1.  Parameters  for  Numerical  Cases 


No. 

Tunnel 

Mo 

Re^/ft 

^OD' 

°R 

Pay 

IWft2 

Tw, 

°R 

hf- 

ft 

a 

e. 

deg 

IMAX 

JMAX 

1 

VKF  B 

7.94 

0.96X1Q6 

93.68 

2.885 

520.0 

0. 035 

a  35 

a 

46 

40 

2 

VKF  B 

7.94 

a96xl06 

93.68 

2.885 

52a  0 

0.035 

a  35 

8.6 

66 

40 

3 

NASA  LRC 
Mach  8 

7.73 

0.48X  106 

100.6 

2. 120 

53a  4 

0.035 

0 

11.1 

60 

40 

D 

VKF  B 

7.93 

0. 97  x  100 

94.43 

3. 140 

525.0 

0,030 

a  25 

10  7 

171 

25 

5 

VKF  B 

7.94 

0. 96  x  1(0 

93.68 

2.885 

520.0 

0.035 

a  35 

a 

181 

35 
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Figure  7.  Flat  plate  velocity  profile. 
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a.  Shock  wave  shape 


008 


b.  Leading  edge  detail 

Figure  9.  Leading  edge  shock  shape  for  a  flat  plate  at  hypersonic 
velocities. 
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Figure  10.  Laminar  hypersonic  boundary— layer/shock— wave  interaction 
using  VKF  Tunnel  B  conditions. 
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An  additional  hypersonic  laminar  boundary-layer/shock-wave 
interaction  using  as  a  basis  the  data  of  Kaufman  and  Johnson  (Ref.  23) 
taken  in  the  Langley  Research  Center  (LRC)  Mach  8  Variable  Density 
Tunnel  was  attempted.  The  results  of  the  computation  as  well  as  the 
relevant  experimental  data  are  presented  in  Fig.  11.  Agreement  is 
satisfactory  particularly  near  the  peak  pressure  location.  The  predicted 
pressure  gradient  within  the  interaction  region  agrees  well  with  the 
experimental  data.  The  major  discrepancy  occurs  near  the  beginning 
of  the  interaction  region  where  the  experimental  pressure  is  about 
50  percent  higher  than  the  computed  pressure.  It  is  tempting  to 
ascribe  the  anomaly  to  viscous-induced  phenomena  from  the  leading 
edge  of  the  plate,  but  the  results  of  Fig.  6  give  some  confidence  in 
the  ability  of  the  code  to  properly  compute  viscous-induced  effects. 

Local  variations  in  wall  temperature  as  well  as  the  overall  model 
temperature  can  affect  the  pressure  in  the  vicinity  of  the  leading 
edge.  Kaufman  and  Johnson  suggest  using  a  uniform  model  temperature 
of  0.5  times  the  total  temperature.  This  suggestion  was  followed 
in  making  the  computations  and  could  have  contributed  to  the  discrepancy 
as  the  exact  model  temperature  distribution  was  not  known.  Six  hours 
of  IBM  370/155  CPU  time  were  required. 

Heat  transfer  is  another  of  the  quantities  commonly  measured  in 
experiments  and  needed  from  calculations.  Thus,  the  capability  of 
the  code  to  predict  the  level  of  heat  transfer  in  a  laminar  boundary- 
layer/shock-wave  interaction  is  also  of  interest.  Figure  12  was 
generated  using  as  a  basis  heat-transfer  data  from  the  VKF  Tunnel  B 
(Ref.  18).  The  agreement  between  the  computed  free-stream  Stanton 
number  (St^)  and  the  measured  free-stream  Stanton  number  is  disap¬ 
pointing.  The  computed  solution  underpredicts  the  maximum  heating 
rate  and  overpredicts  the  minimum  heating  rate.  The  numerical 
solution  indicates  no  separation  to  be  present,  and  an  examination 
of  the  data  indicates  little,  if  any,  separation.  Thus,  the  problem 
appears  to  be  one  of  grid  resolution  as  qualitatively  the  correct 
trends  are  observed  in  the  computer-generated  solution. 
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Figure  11.  Laminar  hypersonic  boundary-layer/si 
using  LRC  Mach  B  tunnel  conditions. 
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Figure  12.  Heat  transfer  for  a  laminar  hypersonic  boundary-layer/shock— wave 
interaction  using  VKF  Tunnel  B  conditions.  . 
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The  cell  Reynolds  number  as  an  index  of  spatial  resolution  has 
previously  been  examined  by  Roache  (Ref.  2),  Cheng  (Ref.  24),  and 
MacCormack  (Ref.  11).  MacCormack,  for  example,  suggests  a  cell 
Reynolds  number  on  the  order  of  two  for  each  coordinate  mesh  spacing 
if  every  term  of  the  Navier-Stokes  equations  is  to  receive  adequate 
support.  For  this  particular  case,  a  cell  Reynolds  number  based  on 
Ay  within  the  viscous-dominated  region  has  a  typical  value  of  6.0;  i.e.. 


ReAy  = 


pv  Ay 
P 


6.0 


(26) 


and  the  cell  Reynolds  number  based  on  Ax  within  the  viscous-dominated 
region  has  a  typical  value  of  1,250;  i.e.. 

Re  a  =  -2^  «  1,250  (  27) 

Thus,  the  y-mesh  should  adequately  support  all  terms,  while  the  x-mesh 
will  not  adequately  support  all  terms.  Obviously,  Ax  cannot  be  decreased 
by  the  several  orders  of  magnitude  needed  for  adequate  support  of  all 
the  terms  in  the  Navier-Stokes  equation  as  CPU  time,  and  core  storage 
would  become  totally  impractical.  MacCormack  (Ref.  11)  suggests  that 
the  inadequately  supported  terms  are  not  necessary  for  boundary-layer 
flow  calculation,  but  anomalies  between  the  experimental  data  and  the 
numerical  solution  show  quantitative  differences,  which  could  be  the 
results  of  inadequate  mesh  resolution  in  the  x-direction. 


Additionally,  at  the  hypersonic  conditions  used  in  this  report,  it 
was  necessary  to  add  the  damping  terms  previously  discussed  to  stabilize 
the  solution.  These  terms,  especially  in  regions  of  large  gradients, 
could  cause  "smearing"  of  the  numerical  solution  as  exemplified  in 
Fig.  12.  It  is  not  possible  to  say  which  of  the  two  effects,  grid 
resolution  or  damping,  caused  the  discrepancies  between  the  experimental 
and  numerical  results.  Approximately  6  hr  of  IBM  370/165  CPU  time 
were  required  for  the  above  solution. 
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The  last  case  to  be  examined  is  based  on  VKF  Tunnel  B  conditions 
and  possesses  a  large  separated  region.  Figure  13  summarizes  the 
tunnel  conditions  used  and  presents  the  results  of  a  time-dependent 
numerical  solution.  This  solution  was  generated  by  computing  a 
hypersonic  flat-plate  solution  and  using  this  solution  as  initial 
conditions  for  the  upstream  portion  of  the  computational  region. 
Agreement  between  computed  and  experimental  data  was  adequate  in  the 
region  about  the  peak  pressure  and  good  for  the  pressure  gradient  near 
the  end  of  the  interaction.  The  middle  portion  of  the  interaction 
exhibits  the  correct  plateau  of  pressure.  The  largest  area  of  dis¬ 
agreement  is  the  beginning  of  the  interaction  where  the  computed 
plateau  and  separation  regions  are  much  shorter  than  those  indicated 
by  experimental  data. 

Several  attempts  were  made  to  rectify  the  situation,  but  none 
were  successful.  A  larger  separation  region  was  used  in  the  initial 
conditions  for  a  subsequent  series  of  runs,  but  as  the  solution  evolved, 
the  separated  region  shrank  to  the  initial  expanse  as  seen  in  Fig.  13. 
Another  attempt  was  made  using  input  profiles  at  0.25  ft  from  the 
leading  edge  instead  of  0.50  ft.  The  final  result  after  some  hours  of 
computer  time  did  not  differ  appreciably  from  that  shown  in  Fig.  13. 


As  in  the  previous  case,  grid  resolution  effects  and  the  damping 
terms  effects  are  cited  as  being  possible  causes  of  the  discrepancies 
between  experimental  data  and  the  results  of  the  numerical  solution. 

The  cell  Reynolds  numbers  and  possess  approximately  the  same 

typical  values  for  this  case  as  in  the  case  previously  examined.  Hence 
the  x-mesh  will  not  adequately  support  every  term  in  the  Navier-Stokes 
equations.  The  damping  terms  were  used  to  stabilize  the  numerical 
solution  generated  for  this  case.  The  gradients  downstream  of  the 
pressure  plateau  region  are  larger  than  the  gradients  upstream  of  the 
plateau  region.  The  values  of  the  damping  terms  were  examined  for  the 
entire  flow  field,  and  in  general,  the  values  in  the  region  downstream 
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of  the  pressure  plateau  were  an  order  of  magnitude  larger  than  the 
values  upstream  of  the  pressure  plateau.  The  magnitude  of  the  damping 
terms  in  the  region  upstream  of  the  pressure  plateau  (in  the  region 
about  the  point  of  separation)  were  larger  than  the  remaining  portion 
of  the  flow  field  except  for  the  previously  delineated  downstream  region. 
Thus,  the  damping  terms  are  of  less  relative  importance  within  the 
separated  recirculation  region,  while  grid  resolution  effects  are  of 
equal  importance.  Close  to  the  point  of  separation,  the  boundary- 
layer  equations  exhibit  singular  behavior,  and  within  the  separated 
recirculation  region,  the  Navier-Stokes  equations  hold.  It  then 
follows  that,  within  this  region,  terms  not  otherwise  of  importance 
may  be  significant.  Lack  of  grid  resolution  in  the  x-direction  appears 
to  be  a  likely  candidate  for  the  anomalous  behavior  of  the  numerical 
solution.  Large  separation  regions  in  laminar  flow  are  difficult  to 
predict  numerically  as  evidenced  by  the  results  of  Skoglund  and  Gay 
(Ref.  6)  and  Hung  and  MacCormack  (Ref.  12).  It  seems  possible  that 
the  difficulties  encountered  with  accurate  prediction  of  large  separation 
region  are  caused  by  inadequate  grid  resolution. 

Although  the  region  of  separation  was  underpredicted  in  extent, 
some  interesting  information  was  obtained  from  the  numerical  solution. 
Figure  14  illustrates  two  velocity  profiles:  one  at  separation  and 
one  in  the  separated  region.  The  profile  at  separation  exhibits  the 
expected  behavior.  The  profile  within  the  separated  region  shows 
appreciable  back  flow,  a  ’’negative1'  wall  shearing  stress,  and  a 
substantial  height  of  the  separation  bubble.  The  enormous  compres¬ 
sion  that  a  boundary  layer  goes  through  within  the  interaction  is 
well  illustrated  by  Fig.  15.  The  streamlines  leaving  the  interaction 
are  compressed  to  a  small  fraction  of  their  height  entering  the 
interaction.  Figure  16  presents  streamline  patterns  for  the  complete 
computational  region.  The  initial  turning  of  the  flow  by  the  incident 
shock  wave  as  well  as  the  location  of  the  strong  shock  wave  downstream 
of  the  separated  region  can  be  seen.  Approximately  16  hr  of  IBM  370/165 
CPU  time  were  required  to  obtain  this  solution. 
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Figure  14.  Velocity  profiles  at  the  point  of  separation  and  within  the  separated  region. 
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Figure  15.  Streamlines  for  hypersonic  laminar  boundary— layer /shock— wave 
interaction. 
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4.0  CONCLUDING  REMARKS  . 

The  method  of  MacCormack  has  been  used  to  solve  numerically  the 
laminar  compressible  time-dependent  Navier-Stokes  equations  for  several 
boundary-layer/shock-wave  interactions  in  the  hypersonic  regime. 

Solutions  for  each  of  the  interactions  were  generated  using  conditions 
corresponding  to  available  experimental  studies.  Comparisons  of  the 
numerical  solutions  with  experimental  data  were  made  to  ascertain  the 
validity  of  the  numerical  method  and  to  identify  regions  of  anomalous 
behavior.  Possible  causes  of  the  anomalies  were  examined. 

The  algorithm  gave  good  results  when  applied  to  laminar  hypersonic 
interactions  that  possessed  small  or  no  separated  regions.  The  pre¬ 
dicted  pressure  distributions  were  in  excellent  agreement  with  experi¬ 
mental  data  for  cases  with  small  or  non-existent  separated  regions. 

The  predicted  wall  heat-transf er  rates  exhibited  the  correct  qualitative 
trends  but  not  the  experimentally  measured  quantitative  values.  Overall, 
the  code's  performance  should  be  considered  as  adequate  for  the  predic¬ 
tion  of  laminar  unseparated  boundary-layer/shock-wave  interactions. 

The  method,  when  applied  to  hypersonic  interactions  having  large 
separated  and  recirculating  regions,  gave  marginal  performance.  The 
region  of  separation  was  underpredicted  by  a  considerable  amount 
although  the  correct  plateau  pressure  and  the  correct  pressure  gradient 
near  the  end  of  the  interaction  were  adequately  predicted.  Considera¬ 
tion  of  the  effects  of  damping  as  well  as  grid  resolution  suggested 
inadequate  mesh  spacing  in  the  x-direction  as  the  cause. 

If  inadequate  mesh  spacing  (and  the  corresponding  lack  of  support  for 
every  term  of  the  Navier-Stokes  equations)  is  a  prime  cause  of  the  cited 
discrepancies  between  the  numerical  results  and  the  experimental  data, 
then  the  needed  reduction  of  several  orders  of  magnitude  in  Ax  would 
increase  CPU  time  and  core  storage  requirements  to  untenable  levels. 
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NOMENCLATURE 

2  2 

C  Constant  volume  specific  heat,  4,290  ft  /sec  -°R 

v 

c  Local  speed  of  sound 

e  Total  internal  energy  per  unit  volume 

F  Column  vector  containing  convection  and  diffusion  fluxes 

in  the  x-direction  defined  by  Eq.  (3) 

G  Column  vector  containing  convection  and  diffusion  fluxes 

in  the  y-direction  defined  by  Eq.  (4) 

h  Height  of  computational  grid,  ft 

hj.  Height  of  fine  mesh  portion  of  computational  grid,  ft 

IMAX  Total  number  of  mesh  points  in  the  x-direction 

JMAX  Total  number  of  mesh  points  in  the  y-direction 

k  Thermal  conductivity 

L  Location  of  shock  impingement  point  on-  flat  plate 

Lx  Operator  denoting  MacCormack  algorithm  x-direction 

contribution  defined  by  Eqs.  (8)  and  (9) 

Ly  Operator  denoting  MacCormack  algorithm  y-direction 

contribution  defined  by  Eqs.  (10)  and  (11) 

M  Fine  grid  passes  for  each  course  grid  pass  defined  by 

Eq.  (25) 
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00 
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Re 

Re 

Re 

St 


Ax 

Ay 


oo 
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AF 

AG 


Free-stream  Mach  number 

Static  pressure 
Reynolds  number 

Cell  Reynolds  number  for  x-mesh  defined  by  Eq.  (27) 

Cell  Reynolds  number  for  y-mesh  defined  by  Eq.  (26) 

Stanton  number  based  on  free-stream  conditions  and 
adiabatic  wall  temperature 

Static  temperature 

Time 

Column  vector  of  conserved  quantities  per  unit  volume  defined 
by  Eq.  (2) 

Component  of  velocity  in  the  x-direction 
Component  of  velocity  in  the  y-direction 
Coordinate  along  plate  surface 
Coordinate  normal  to  plate  surface 

Damping  term  in  Lx  operator  used  to  stabilize  calculation 
defined  in  Eq.  (19) 

Damping  term  in  Ly  operator  used  to  stabilize  calculation 
defined  in  Eq.  (18) 
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At 


Time  step  used  in  temporal  advancement 


Ax 


Mesh  spacing  in  x-direction 


Ay 


Mesh  spacing  in  y-directlon 


Ratio  of  specific  heats,  1.4 


Coefficient  of  bulk  viscosity,  taken  as  -2/3  u 


Molecular  viscosity  coefficient 


Mass  density 


ayLtay  Normal  stress  fluxes  in  Navier-Stokes  equations  defined  by 
Eqs*  (5)  and  (7) 


t  ,T  Shear  stress  fluxes  in  Navier-Stokes  equations  defined  by 
xy  yx 

Eq.  (6) 


Subscripts 


Upstream  of  incident  shock 


Downstream  of  incident  shock  and  upstream  of  reflected  shock 


Downstream  of  reflected  shock 


Denotes  coarse  grid 


Denotes  fine  grid 


General  grid  point  in  x-direction 
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j  General  grid  point  in  y-direction 

w  Evaluated  at  wall 


x  x-direction 


y  y-direction 

o  Stagnation  or  total 

00  Free  stream 


Superscript 

n  Time  level 


Evaluated  during  corrector  pass  using  predictor  values 
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